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1.0  INTRODUCTION 


The  MACH2  (Ref.  1)  two-dimensional  magnetohydrodynamics  (MHD)  computer  code  is  an 
important  tool  for  the  research  conducted  in  the  High  Energy  Plasma  Branch  of  the  Phillips 
Laboratory.*  Like  all  non-Lagrangian  fluid  simulation  codes,  its  accuracy  is  dependent  on  the 
treatment  of  the  advection  terms  of  the  difference  equations.  These  terms  carry  physical 
properties  such  as  density,  internal  energy  and  momentum  with  the  fluid  as  it  moves  across  the 
computational  grid.  MACH2  may  be  run  in  a  purely  Lagrangian  mode,  where  the  grid  vertices 
move  with  the  fluid,  but  this  is  rarely  a  good  approach.  Lagrangian  computational  grids  usually 
become  entangled  when  the  fluid  has  variations  in  more  than  one  dimension,  and  this  causes  the 
code  to  crash.  MACH2  has  an  adaptive  grid  generator  (Refs.  2  and  3)  which  can  be  used  to 
create  an  almost  Lagrangian  grid  that  remains  smooth.  However,  any  grid  that  is  not  purely 
Lagrangian  will  depend  on  the  advection  algorithm.  Further,  using  a  purely  Eulerian  grid,  whose 
vertices  remain  at  fixed  positions,  is  often  the  easiest  approach  for  the  initial  simulations  in  a 
study,  and  it  may  be  the  only  practical  rqiproach  for  a  veiy  complicated  geometry.  This  report 
discusses  an  in-house  project  to  replace  the  first-order  accurate  Godunov  advection  algorithm 
(Ref.  4)  in  MACH2  with  a  more  accurate  van  Leer  algorithm  (Ref.  5). 

Fluid  simulation  codes  solve  a  set  of  coupled  partial  differential  equations  that  are  conservation 
statements  for  mass,  momentum  and  energy.  The  solution  procedure  for  these  equations  in 
MACH2  is  an  Arbitrary-Lagrangian-Eulerian  (ALE)  algorithm  (Ref.  6).  During  each  time  step, 
it  solves  for  changes  in  the  physical  quantities  on  a  Lagrangian  grid,  creates  a  new  grid  with  the 
adaptive  grid  generator,  and  then  applies  the  advection  terms  separately  to  meq)  the  Lagrangian 
results  onto  the  new  grid.  Separating  the  advection  terms  is  not  a  new  idea— see,  for  example. 
References  6  and  7,  but  the  adjqitive  grid  generator  gives  MACH2  much  more  flexibility  than 
most  codes. 


♦This  code  was  developed  under  contract  for  the  Air  Force  Weapons  Laboratory— currently 
Phillips  Laboratory— by  Mission  Research  Corporation.  It  features  a  computational  block 
structure  that  allows  a  user  to  model  complicated  geometries  with  ease.  It  has  been  applied  to  a 
large  number  of  Air  Force  projects  including  plasma  flow  switches,  cylindrical  implosions, 
plasma  guns,  and  plasma  toroid  experiments,  and  it  has  been  ^qiplied  to  many  other  Department 
of  Defense  projects. 
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This  report  concentrates  on  the  advection  terms  for  mass,  internal  energy  and 
momentum.  Because  MACH2  treats  conducting  fluids,  it  solves  Faraday’s  law  for  the 
evolution  of  the  magnetic  field,  and  the  appropriate  form  of  this  law  has  an  advection 
term.  However,  this  term  requires  special  attention  and  will  be  addressed  in  a  separate 
report.  Section  2.0  of  this  report  gives  the  analytical  and  numerical  formulations  of  the 
advection  terms.  Section  3.0  explains  how  the  numerical  formulation  is  implemented  in 
MACH2.  Section  4.0  discusses  test  problems  that  illustrate  the  advantages  of  the 
algorithm.  The  final  section  concludes  the  report  with  a  discussion  of  how  the  algorithm 
has  influenced  simulations  of  real  problems. 
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2.0  MATHEMATICAL  DESCRIPTION 


2.1  ANALYTIC  FORMULATION 


The  Transport  Theorem  can  be  used  to  find  the  analytic  relation  between  the  time-derivative  of 
Lagrangian  volume-integrals  and  the  time-derivative  of  other  volume-integrals.  Marsden  and 
Tromba  (Ref.  8)  state  the  scalar  and  vector  forms  of  the  theorem.  Their  definitions  may  be 
expanded  in  two  ways  without  changing  the  results.  They  define  a  time-independent  velocity 
vector  field  which  describes  the  motion  of  fluid  elements  and  integration  domains  that  are 
composed  of  these  elements.  First,  the  velocity  vector  field  may  be  time-dependent  because  it  is 
not  differentiated  with  respect  to  time  in  the  proof.  Second,  the  velocity  vector  field  may 
describe  an  imaginary  fluid— one  that  is  not  moving  widi  the  physical  fluid,  because  no  physical 
laws  are  applied.  For  a  pertinent  example  of  an  imaginary  fluid,  consider  the  time-dependent 
vertex  positions  generated  by  an  adaptive  mesh  algoridim  to  be  markers  on  an  imaginary  fluid. 
With  these  definitions,  the  scalar  and  vector  forms  of  the  theorem  are,  respectively 
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where  f  is  a  scalar  function  of  position,  x,  and  time,  t,  and  W  and  G  are  vector  functions  of 
position  and  time.  The  functions  f  and  G  are  arbitrary,  but  W  is  the  velocity  vector  field  of  the 
real  or  imaginary  fluid  that  carries  the  moving  region  C>.  The  second  and  third  terms  in  the 
integrand  on  the  right  side  of  Equations  1  and  2  may  be  combined  into  the  divergence  of  a 
vector,  f  W,  for  the  former  and  the  divergence  of  a  tensor,  W  f  G,  for  the  latter.  The 
Divergence  Theorem  may  be  applied  to  each  equation  with  the  following  results 
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where  n  is  the  outward  normal  on  the  closed  surface  dC>. 


(3) 


(4) 


If  W  is  the  physical  fluid  velocity,  V,  then  each  of  the  above  expressions  relates  the  time 
derivative  of  a  moving  Lagrangian  volume  to  the  time-derivative  of  a  corresponding  stationary 
volume  plus  a  surface  flux  term.  Consider  also  the  same  two  equations  with  U,  an  imaginary 
fluid  velocity,  substituted  for  W.  The  resulting  expressions  are  related  to  the  same  stationary 
volume  integrals  as  the  Lagrangian  equations,  provided  that  the  physical  fluid  region  and  the 
imaginary  fluid  region  correspond  at  t  =  t.  Subtracting  the  scalar  relation  (Eq.  3)  for  the 
Lagrangian  volume  from  the  same  relation  for  the  volume  of  imaginary  fluid  results  in  the 
following 
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where  H  is  the  Lagrangian  volume  moving  with  the  physical  fluid,  and  Y  is  the  volume  moving 
with  the  imaginary  fluid.  The  vector  relation  obtained  from  Equation  4  is 
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2.2  NUMERICAL  FORMULATION 

The  difference  forms  of  Equations  S  and  6  are  the  equations  solved  in  the  advection  step  in 
MACH2  using  U  as  the  velocity  field  of  the  grid  as  suggested  in  the  example  above.  The 
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difference  form  of  the  time-derivative  terms  are  simple.  Each  is  the  difference  between  a  new 
volume  integral  and  an  old  volume  integral,  divided  by  the  time  increment,  and  a  volume  integral 
is  the  average  of  a  quantity  in  a  cell,  multiplied  by  the  volume  of  the  cell.  Because  the 
Lagrangian  cells  and  the  mesh  cells  are  equivalent  prior  to  the  Lagrangian  step,  these  "old" 
volume  integrals  may  be  eliminated  from  the  difference  equations.  The  Lagrangian  step  of  the 
momentum  algorithm  produces  the  "new"  averages  in  the  Lagrangian  cells,  so  the  new 
Lagrangian  volume  integrals  are  known  at  the  start  of  the  advection  step.  After  multiplying  by 
the  time  increment,  At,  the  following  difference  form  of  Equation  5  results. 

where  f  is  the  new  average  value  of  f  in  the  cell  and  t)  is  the  new  ceU  volume  for  either  the 
Lagrangian  ceU,  where  the  index  ^>  =  fl,  or  the  grid  cell,  where  O  =  ^.  The  summation  on  the 
right  side  represents  the  flux  of  f  over  the  cell  surface  for  one  advection  step.  The  sides  are 
indexed  by  s,  and  each  has  a  surface  area,  a.  The  difference  form  of  Equation  6  is 

(f  G  t)  )  =  (f  G  \)  )  +  Atl<a  n  (  U-  V)f  G>  (8) 

where  G^  is  the  new  cell  average  of  G.  After  solving  for  the  left  side  of  Equation  7  or  8,  one 
may  divide  by  the  new  grid  cell  volume  to  obtain  the  desired  new  grid  cell  average.  As  van  Leer 
points  out,  these  equations  are  exact  (Ref.  5),  and  the  accuracy  of  an  algorithm  depends  on  the 
average  flux  terms  inside  the  summation.  Note  that  when  the  grid  moves  with  the  fluid,  U  =  V, 
and  the  Lagrangian  values  remain  unchanged.  The  grid  is  stationary  when  U  =  0,  and  in  this 
case,  the  two  equations  form  an  Eulerian  representation. 

To  understand  the  van  Leer  approach,  it  is  easiest  to  start  with  the  Godunov  or  ’donor  cell’ 
method.  For  the  latter,  the  average  flux  terms  in  Equations  7  and  8  are  rather  single.  Many 
fluid  codes  have  velocities  centered  between  cells,  so  Newton’s  Second  Law  may  be  written  in  a 
centered  difference  form  with  pressures  being  cell-centered  (along  with  mass  and  internal 
energy).  Thus,  the  velocities  and  areas  in  the  flux  terms  are  located  at  the  interfaces.  To 
construct  the  entire  flux  term,  Godunov  considers  the  cell-centered  quantities  to  be  constant 
within  each  cell.  In  one-dimension  the  representation  is  a  set  of  slabs,  and  one  example  is 
illustrated  in  Figure  1 .  The  flux  at  each  interface  becomes  the  interface  velocity  multiplied  by 
the  magnitude  of  the  slab  from  which  the  flow  comes,  hence  the  name  ’donor  cell.  ’  The  method 
is  explicit  because  the  velocities  and  slab  magnitudes  used  in  the  flux  terms  are  those  at  the 
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Figure  1.  Godunov  representation  of  a  cell-centered  quantity  on  a  on-dimensional  grid. 

beginning  of  the  advection  step.  It  is  also  conservative  because  the  flux  that  comes  out  of  the 
donor  cell  goes  into  the  adjacent  cell. 

When  describing  the  explicit  nature  of  a  method,  one  is  tempted  to  refer  to  the  beginning  of  the 
time  step,  but  this  would  be  misleading.  Simulation  codes  often  treat  physical  processes 
independently  during  each  time  step.  This  suggests  a  linearization  in  time,  which  is  reasonable 
as  long  as  the  time  step  is  small.  In  MACH2  the  changes  from  each  phyisical  process  are  added 
sequentially.  Furthennore,  for  any  ALE  code,  the  advection  step  is  a  mapping  from  the 
Lagrangian  grid,  so  the  Lagrangian  cell-centered  values  are  the  explicit  values  during  the 
advection  step. 

To  obtain  more  accuracy  than  the  Godunov  method,  the  van  Leer  technique  replaces  the  slab 
approximation  of  the  distribution  with  a  better  approximation  (Ref.  5).  It  uses  derivatives  to 
make  the  distribution  a  set  of  trapezoids  instead  of  slabs  (Fig.  2).  This  representation  is 
piecewise  continuous,  like  the  slab  representation,  and  the  correct  cell  averages  are  preserved. 


Figure  2.  Van  Leer  representation  of  a  cell-centered  quantity  on  a  one-dimensional  grid. 
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Van  Leer  proposes  several  possible  formulations  for  the  derivatives.  In  the  simplest  scheme,  the 
upwind  cell-centered  quantity  and  its  corresponding  centered  difference  are  used  to  construct  the 
flux  terms.  Any  subsequent  mention  of  the  "centered-difference  scheme"  refers  to  this  approach. 
The  scalar  flux  term  from  Equation  7  with  this  scheme  is 

(a  n  (  U-  V)f)  =aC 

i  8 

where  f  is  the  upwind  cell-centered  quantity,  Af^  is  its  centered  difference  with  respect  to  the 
advection  direction.  Ax  is  the  cell  dimension  in  the  advection  direction,  and 
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The  sign  of  the  Af  term  in  Equation  9  is  positive  when  the  s-interface  is  the  upper  bound  of  the 
upwind  cell  and  negative  when  it  is  the  lower  bound.  When  the  time  step  is  limited  by  the 
Courant-Friedrichs-Lewy  (CFL)  condition,  |  ^  (At/Ax)  |  ^  1,  the  term  in  the  parentheses  on  the 
right  side  is  bounded  by  zero  and  one,  so  this  scheme  limits  to  the  donor  cell  method  as 
I  ^^(At/Ax)  I  approaches  unity.  Tlie  vector  equation  corresponding  to  Equation  9  has  f  G  in 
place  of  f.  For  multiple  dimensions,  there  will  be  multiple  centered  differences  for  each  quantity. 

Van  Leer  derives  the  an^lification  factor  for  this  scheme  with  uniform  grid  spacing  and 
velocities. 
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where  O  =  ^(At/Ax),  and  the  index  is  omitted  to  indicate  the  uniformity  (Ref.  S).  The  angle 
a  =  27cAx/1,  and  1  is  the  length  of  a  wave  moving  across  the  grid.  The  dissipation  error  per  time 
step  is  one  minus  the  magnitude  of  the  amplification  factor,  and  it  has  a  maximum  at  o  =  1/2 
(Ref.  5).  Thus,  for  0=1/2  and  a  =  n/2,  the  centered-difference  scheme  has  a  dissipation  error  of 
0.12  per  time  step.  The  amplification  factor  for  the  Godunov  method  is 
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(11) 


=  1  -  a(  1  -  cosa )  -  i  o  sina 

Its  dissipation  error  for  a  =  1/2  and  a  =  nfl  is  O.S  per  time  step,  so  even  the  simplest  van  Leer 
scheme  is  a  considerable  improvement  for  intermediate  length  waves.  The  dispersion  error  is 
measured  by  the  ratio  of  the  numerical  advection  speed  to  the  ‘me  advection  speed, 

(0  =  arg(g)  /  (-oa).  Polar  plots  of  I  g  I  at  o  =  1/2  and  (O  in  the  limit  of  vanishing  a  for  both  the 
centered-difference  scheme  and  the  Godunov  method  are  shown  in  Figure  3.  A  sufficient 
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(a)  Magnitude  of  ampliflcation  factors  at  a  =  1/2. 


(b)  Ratio  of  numerical  advection  speed  to  true  advection  speed  in  the  limit  of  vanishing  a. 


Figure  3.  Performance  of  the  advection  methods.  The  solid  line  illustrates  the  centered 
difference  scheme,  and  the  dashed  line  illustrates  the  Godunov  method.  The 
angle  from  the  positive  horizontal  axis  is  a  =  2»rAx/l. 


8 


condition  for  stability  is  that  I  g  I  ^  1 ,  for  if  this  were  not  satisfied,  rq)eated  qiplications  of  the 
advection  step  would  force  wave  amplitudes  to  grow  geometrically.  Both  the  Godunov  and 
centered-difference  scheme  meet  this  criterion  when  the  CFL  condition  is  satisfied.  Although  it 
is  not  accurate  to  make  generalizations  about  the  stability  of  difierence  equations,  the  upwind 
nature  of  an  advection  scheme  tends  to  add  stability,  and  both  of  these  schemes  have  an  upwind 
nature. 

Besides  accuracy  and  stability,  the  monotonicity  of  a  distribution  should  be  respected  by  the 
advection  algorithm.  To  quote  van  Leer  (Ref.  5),  "The  monotonicity  condition  says  that,  when  a 
monotonic  initial  value  distribution  is  numerically  convected,  the  resulting  distribution  must  be 
monotonic  again."  Hiis  is  enforced  by  placing  limits  on  the  centered  difference,  Af^,  in 
Equation  9.  The  limits  prevent  the  linear  distribution  of  a  quantity  within  a  cell  from  exceeding 
the  cell-centered  average  of  that  quantity  in  the  adjacent  cells.  Also,  if  a  cell-centered  average  is 
not  between  those  of  the  adjacent  cells,  the  slab  representation  is  used.  This  prevents  the 
development  of  new  extrema.  The  representation  in  Figure  2  is  properly  limited. 

The  centered-difference  scheme  requires  only  a  small  amount  of  additional  computation  time 
over  the  Godunov  method,  most  of  which  is  spent  on  finding  the  centered  differences.  These 
differences  are  calculated  during  the  advection  algorithm  and  do  not  require  permanent  storage. 
The  other  second-order  schemes  proposed  by  van  Leer  are  based  on  derivatives  that  are 
computed  separately  from  the  cell-centered  quantities.  These  derivatives  require  separate 
storage,  and  they  must  be  updated  during  all  of  the  other  algorithms  that  change  the 
corresponding  cell-centered  quantities.  The  accunu;y  analysis  in  Reference  S  shows  that  some 
of  the  more  complicated  schemes  can  track  waves  as  riiort  as  two  cell  lengths  widi  very  little 
dissipation  and  dispersion,  whereas  the  scheme  with  the  centered  differences  loses  accuracy  for 
wavelengths  less  than  four  cell  lengths.  However,  for  enhancing  an  existing  fluid  code,  it  is  far 
easier  to  add  the  centered-difference  scheme  to  the  advection  algorithm  than  it  is  to  rewrite  the 
entire  code  to  track  derivatives.  Therefore,  the  centered-difference  scheme  has  been  added  to 
MACH2. 
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3.0  IMPLEMENTATION 


3.1  PRELIMINARIES 

For  a  simple  one-dimensional  advection  problem  with  uniform  velocities  and  surface  areas 
between  cells,  the  terms  of  Equation  9  have  been  sufficiently  deffiied.  For  anything  more 
complicated,  they  are  ambiguous.  The  resolution  of  the  ambiguities  is  a  code-dependent  issue. 
This  section  wiU  address  this  issue  for  version  v9101  of  MACH2,  and  it  will  provide  a  guided 
tour  through  the  subroutines  for  those  who  use  and  modify  the  code. 

The  subroutine  ARUN  contains  the  main  loop  of  MACH2  that  calls  separate  routines  for  each  of 
the  physical  processes.  The  subroutine  HYDRO,  which  is  called  from  ARUN,  calculates  the 
Lagrangian  stage  of  the  ALE  algorithm.  Following  HYDRO,  ARUN  calls  REMESH.  The  first 
part  of  this  subroutine  calls  the  ad2q)tive  mesh  generator,  and  the  second  part  calls  TRNSPT,  the 
advection  algorithm,  to  complete  the  ALE  algorithm. 

When  MACH2  is  used  for  planar  geometries,  the  x-y  plane  is  the  computational  plane,  and  no 
variations  are  allowed  in  the  perpendicular  direction.  For  axisymmetric  geometries,  the  r-z  plane 
forms  the  computational  domain.  For  either  case,  TRNSPT  calculates  the  advection  that  results 
from  velocity  components  in  the  computational  plane.  Other  advection  terms  that  result  from  the 
0-component  of  velocity  are  also  nonzero  in  axisymmetric  geometries.  These  terms  are  treated 
at  the  end  of  HYDRO  and  are  not  considered  in  TRNSPT. 

To  simulate  complicated  geometries  with  MACH2,  the  spatial  domain  is  decomposed  into 
four-sided  blocks— the  reader  is  encouraged  to  see  Reference  1  for  more  information  on  allowed 
domains.  The  blocks  are  divided  into  quadrilateral  cells  which  form  the  computational  grid. 
Cells  have  a  horizontal  index,  i,  and  a  vertical  index,  j.  Each  physical  algorithm  solves  or 
iterates  its  process  on  one  block  at  a  time,  and  boundary  conditions  couple  adjacent  blocks.  The 
spatially-dependent  physical  quantities  are  stored  in  "pointered"  memory  location;  i.  e.,  the 
POINTER  extension  of  standard  FORTRAN  is  used  to  set  the  two-dimensional  arrays  for  the 
physical  quantities  to  the  {q}propriate  set  of  memory  locations  for  each  block.  This  conserves 
memory  because  the  dimension  of  pointered  arrays  can  vary  from  block  to  block.  However,  the 
pointers  and  dimensions  must  be  set  before  the  arrays  can  be  correctly  accessed.  Therefore,  all 
of  the  physical  algorithms  will  contain  "do-loops"  over  the  blocks,  and  the  first  call  is  always  to 
the  subroutine  SETBLK  which  sets  the  pointers  and  dimensions.  These  loops  will  be  mentioned 
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frequently  in  the  description  of  the  advection  algorithm  below. 


3.2  CELL-CENTERED  OUANTmES 

The  advection  of  the  cell-centered  quantities  is  not  complicated.  The  frrst  block  loop  of 
TRNSPT  calls  the  subroutine  TRNSINIT— see  the  Appendix  for  a  listing  of  the  nonmagnetic 
advection  subroutines.  The  first  loop  of  TRNSINIT  creates  the  Lagrangian-cell  volume 
integrals,  each  of  which  forms  the  first  term  on  the  right  side  of  Equation  7  or  8  for  the  volume 
of  one  cell.  For  example,  when  mass  is  advected,  this  integral  is  simply  the  total  mass  of  the 
Lagrangian  cell.  This  is  also  the  same  as  the  mass  of  the  cell  prior  to  the  Lagrangian  phase, 
which  is  the  old  mass  density,  stored  in  the  "ro"  array,  mult^lied  by  the  old  cell  volume, 
"oldvol."  For  the  internal  energy,  the  integral  is  the  total  internal  energy  in  the  cell.  This  is  the 
Lagrangian  cell  mass  just  computed,  "mp,"  multiplied  by  the  specific  internal  energy  after  the 
Lagrangian  stage,  "sel."  This  loop  also  initializes  die  Lagrangian  densities  where  necessary.  For 
internal  energy,  this  quantity  is  the  internal  energy  per  unit  volume.  It  is  found  by  multiplying 
the  Lagrangian  mass  density,  "rol,"  by  "sel,"  and  is  stored  back  in  "sel"  array. 

This  TRNSINIT  loop  also  defines  the  Lagrangian  cell  volumes,  the  "lagvol"  array,  and  the 
volumes  exchanged  between  these  cells  during  the  m^ing  to  the  new  grid.  The  exchange 
volumes  are  determined  with  the  cross  product  of  two  vectors,  the  grid  velocity  relative  to  die 
fluid  and  the  displacement  vector  from  one  vertex  of  the  new-grid  cell  to  its  next  vertex.  The 
relative  velocities  define  the  U  -  V  vector  in  Equation  9  and  are  the  vertex-centered 
("urr',"vrl")  array  pair.  Figure  4  illustrates  these  vectors  for  the  bottom  exchange  volume, 
"dxbott."  The  dimensions  on  the  cross  products  are  area  per  unit  time.  They  are  multq>lied  by 
the  time  increment  "dt"  and  an  iq>propriate  perpendicular  dimension  to  form  a  volume.  For 
planar  geometries,  the  radius  array,  "r,"  is  set  to  unity,  so  the  volume  is  per  unit  depth 
perpendicular  to  the  computational  plane.  For  axisymmetric  geometries,  the  radius  is  factored 
into  the  relative  velocities  to  create  a  volume  per  unit  angle— see  the  TRNSINIT  listing  in  the 
Appendix  for  the  formulation.  These  exchange  volumes  take  the  place  of  the  a^  factor  on  the 
right  side  of  Equation  9. 

The  "200"  loop  of  TRNSINIT  separates  the  "con2"  fraction  of  marker  material  to  advect  its  mass 
separately  from  the  rest  of  the  mass.  The  "300"  and  "400"  loops  compare  the  size  of  the 
exchange  volumes  to  the  Lagrangian  and  new-grid  volumes  and  saves  die  largest  ratio  for  the 
time  step  control. 
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Figure  4.  Illustration  of  the  "dxbott"  exchange  volume  between  a  Lagrangian  cell  and 
a  new  grid  cell. 


Once  the  TRNSINIT  loop  of  TRNSPT  is  complete,  the  cell-centered  quantities  are  separately 
passed  into  the  subroutine  TRNSLP.  The  first  block  loop  in  TRNSLP  calls  TRNSGR  which 
finds  the  centered  differences  of  the  quantity  passed  into  the  routine.  Note  that  these  differences 
are  differences  of  the  Lagrangian  densities.  The  "200"  loop  of  TRNSGR  finds  the  centered 
difference  in  the  j-index  direction,  and  the  "300"  loop  finds  the  difference  in  the  i-index 
direction.  Each  are  limited  to  twice  the  corresponding  backward  and  forward  differences  for 
monotonicity,  and  if  the  signs  of  those  differences  do  not  agree,  the  centered  difference  is 
reduced  to  zero.  This  implementation  is  die  monotonicity  algorithm  of  Equation  66  in 
Reference  5.  The  more  conservative  algorithm  of  Equation  67  in  Reference  5  is  in  the  current 
version  of  MACH2,  v9101. 

The  second  block  loop  of  TRNSLP  calls  TRNSDQBC,  which  communicates  the  differences 
along  the  boundaries  of  adjacent  blocks,  and  subsequendy  calls  TRNSADV,  which  perfotms  the 
advection.  The  "  100"  loop  of  TRNSADV  advects  the  quantity  in  the  j-index  direction,  and  the 
"200"  loop  advects  it  in  the  i-index  direction.  They  create  the  flux  of  the  quantity  on  the  bottom 
side  and  left  side  of  each  cell,  respectively.  This  is  the  application  of  Equation  9  for  scalar 
quantities.  For  vector  quantities,  each  component  is  separately  passed  into  TRNSLP.  The  ratio 
of  the  volume  flux  to  the  donor-cell  Lagrangian  volume  in  TRNSADV  replaces  the  ^  (At/Ax)  in 
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Equation  9.  Both  ratios  represent  the  fraction  of  the  cell  advected,  but  the  volume  ratio 
automatically  accounts  for  arbitrary  cell  shapes.  After  finding  the  flux,  these  loops  remove  it 
from  the  donor-cell’s  Lagrangian  volume  integral  and  add  it  to  the  adjacent  cell’s  integral.  When 
the  TRNSADV  loop  is  complete,  the  integral  array  holds  the  new-grid  integrals,  each  of  which  is 
the  left  side  of  Equation  7  or  8  for  the  new-grid  cell. 

3.3  MOMENTUM 

The  advection  of  momentum  is  more  difficult  than  the  advection  of  cell-centered  quantities.  The 
velocities  are  centered  at  the  grid  vertices  and  not  the  cell  centers,  and  to  use  the  cell-centered 
scheme,  one  must  first  create  cell-centered  momenta.  Margolin  and  Beason  have  suggested 
creating  cell-centered  quantities  that  are  the  average  and  derivatives  of  die  surrounding  vertex 
quantities  (Ref.  9).  The  current  implementation  in  MACH2  is  similar  to  this  q>proach.  For  each 
velocity  component,  it  creates  four  cell-centered  momentum  densities  which  are  the  products  of 
each  of  the  four  vertex  velocities  and  the  cell-centered  mass  density.  Following  the  advection, 
the  four  resulting  cell  integrals  are  distributed  into  vertex  momenta.  To  create  the  new  vertex 
velocity,  each  vertex  momentum  is  divided  by  the  vertex  mass,  which  is  the  average  of  the  four 
cell  masses  surrounding  the  vertex.  This  is  similar  to  the  scheme  in  Reference  9,  for  if  one  used 
the  average  and  three  possible  differences  instead  of  the  average  and  three  possible  derivatives, 
the  advection  is  algebraically  equivalent  to  the  MACH2  scheme. 

Unlike  the  cell-centered  quantities,  TRNSLP  is  not  called  directly  from  TRNSPT.  Instead, 
TRNSPT  calls  the  subroutine  TRNSMM.  This  subroutine  first  calls  TRNSMMIN  to  create  die 
four  momentum  densities  and  Lagrangian  momentum  integrals  for  each  conqionent.  TRNSMM 
then  calls  TRNSLP  for  each  of  the  four.  Before  returning  to  TRNSPT,  the  subroutine 
TRNSMMF  is  called  to  distribute  the  new-grid  momentum  integrals  among  the  vertices.  The 
new  velocities  are  calculated  after  TRNSPT  during  the  REMESH  call  to  RMSHVEL. 

The  monotonicity  of  the  momentum  advection  also  needs  special  attention.  Although  the 
monotonicity  algorithm  in  TRNSGR  will  not  create  new  extrema  in  the  resulting  momentum 
distribution,  it  does  not  guarantee  the  same  for  the  resulting  velocity  distribution.  Consider  a 
situation  where  the  mass  density  distribution  is  monotonic,  and  the  gradient  is  in  the  direction  of 
a  uniform  flow  in  the  computational  plane.  In  addition,  there  is  momentum  density 
perpendicular  to  the  conqiutational  plane,  and  its  distribution  is  not  monotonic.  The  centered 
difference  will  be  used  for  the  advection  of  the  mass  and  the  fraction  of  mass  removed  from  the 
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donor  cell  will  be  larger  than  the  fraction  of  volume  removed.  However,  the  fraction  of 
perpendicular  momentum  removed  from  the  donor  cell  will  be  equivalent  to  the  volume  fraction, 
and  the  resulting  velocity  distribution  will  show  a  new  maxima.  The  results  of  a  simulation  with 
these  conditions  is  presented  in  Section  4.0.  It  has  been  found  that  when  new  maxima  develop  in 
the  velocity  components  that  are  in  the  computational  plane,  a  numerical  instability  may  result. 

One  prescription  to  avoid  the  instability  is  to  discard  the  difference  of  the  momentum  density  and 
use  the  difference  of  the  mass  density,  multiplied  by  the  vertex  velocity,  instead.  With  this 
prescription,  the  amount  of  momentum  advected  is  proportional  to  the  amount  of  mass  advected. 
Mathematically,  this  uses  the  product  rule  on  the  derivative  of  the  momentum, 

d  (pv)  3  V  3  p 

- =p -  +v -  (12) 

3x  3x  3x 

where  p  is  the  mass  density  and  v  is  a  velocity  component,  and  throws  away  the  first  term  on  the 
right  before  convening  to  a  difrerence  form.  This  product  is  formed  during  the  call  to 
TRNSMMGR  in  the  second  block  loop  of  TRNSLP. 

Another  possible  prescription  to  avoid  the  creation  of  new  velocity  extrema  is  to  difference  both 
terms  on  the  right  of  Equation  12  and  apply  the  monotonicity  algorithm  to  each  term  separately. 
This  scheme  has  also  been  attempted  with  MACH2,  and  although  it  advects  velocity  gradients 
with  less  diffusion,  it  seems  to  be  noisy.  One-dimensional  advection  problems  develop  enough 
noise  to  have  noticeable  two-dimensional  variations. 
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4.0  RESULTS 


There  are  two  one-dimensional  test  problems  that  are  often  used  to  evaluate  the  performance  of 
an  advection  algorithm.  The  first  is  the  uniform  advection  of  a  square  pulse  of  some  quantity, 
and  the  second  is  the  shock  tube.  The  Fourier  Transform  of  a  square  pulse  is  an  oscillating, 
continuous  function,  and  the  magnitude  of  the  oscillations  is  inversely  proportional  to  the  wave 
number  and  does  not  diminish  abruptly.  For  the  test  problem,  the  initial  pulse  is  formed  with  an 
integral  number  of  cells.  Thus,  when  the  initial  pulse  has  relatively  few  cells,  the  problem  will 
exercise  an  algorithm’s  ability  to  advect  high  frequency  waves  in  a  manner  that  is  relevant  to 
simulations  of  actual  experiments  and  physical  phenomena. 

The  pulse  test  problem  described  here  is  initiated  in  the  following  manner.  The  domain  is  a  long 
rectangular  chamber,  0.4  m  in  width  and  12.8  m  in  length,  which  is  divided  into  512  square  cells. 
The  fluid  is  given  a  uniform  velocity  component  of  1  m/s  in  the  long  dimension.  The  pulse  is 
positioned  from  0.4  m  to  0.8  m  from  the  left  side  of  the  chamber,  so  it  is  initially  four  cells  long. 
It  is  composed  of  mass  that  is  10  kg/m^  in  density,  which  is  a  factor  of  10  larger  than  the  density 
in  the  rest  of  the  chamber.  It  is  given  a  temperature  of  10’*^  eV,  also  a  factor  of  10  above  that 
outside  the  pulse.  The  extreme  temperatures  are  chosen  to  make  the  sound  speed  very  small 
compared  with  the  advection  velocity.  Finally,  the  time  step  is  limited  so  that  a  ^  1/2. 

Figures  5  and  6  show  the  mass  density  and  temperature  distributions,  respectively,  for  the 
centered-difference  scheme  with  monotonicity  and  the  Godunov  method,  after  the  pulses  travel 
across  100  cells.  Although  both  algorithms  diffuse  the  small  pulse,  the  peak  mass  with  the 
centered-difference  scheme  is  twice  that  obtained  with  the  Godunov  method,  and  the  pulse  width 
is  much  less.  The  temperature  pulses  have  a  different  shtq>e  from  the  density  pulses  and 
maintain  relatively  larger  peaks.  This  occurs  because  the  temperature  is  the  quotient  of  two 
advected  quantities,  internal  energy  and  mass.  In  this  case,  the  initial  internal  energy  pulse  is 
two  orders  of  magnitude  larger  than  the  background.  Note  that  the  algorithm  will  not  create  new 
extrema  in  the  temperature  distribution  because  the  monotonicity  algorithm  prevents  the  creation 
of  new  extrema  in  both  the  internal  energy  and  mass  distributions. 

A  variation  of  the  pulse  problem  can  illustrate  the  difficulty  with  momentum.  Consider,  again,  a 
rectangular  chamber,  with  a  mass  density  of  1  kg/m^  in  half  of  the  chamber  and  10'^  kg/m^  in  the 
other  half.  The  velocity  in  the  computational  plane  is  a  uniform  1  m/s  towards  the  side  with  the 
greater  density,  and  the  less  dense  side  has  a  perpendicular  velocity  component  of  1  m/s.  Thus, 
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(a)  Godunov  method. 


(b)  Centeied-difference  scheme. 


Figure  S.  Mass  density  distributions  for  the  square  pulse  advection. 


(a)  Godunov  method. 


(b)  Centeied-difference  scheme. 


Figure  6.  Temperature  distributions  for  the  square  pulse  advection. 
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the  density  and  perpendicular  velocity  distributions  are  both  step  functions,  but  the  change  of  one 
is  opposite  the  other.  This  is  the  problem  described  in  Subsection  3.3.  Figure  7  shows  two 
distributions  of  perpendicular  velocity— the  distribution  on  the  left  results  when  the  momentum 
difference  is  created  from  the  product  of  the  cell-centered  mass  and  vertex-centered  velocity,  and 
the  distribution  on  the  right  results  when  one  uses  the  density  difrerence  multiplied  by  the 
velocity.  The  new  maximum  in  the  former  is  obvious.  Note  that  the  distribution  of 
vertex-centered  perpendicular  momentum,  defined  by  the  product  of  the  vertex-centered  velocity 
component  and  an  average  mass  of  the  adjacent  cells,  has  a  maximum  at  the  begirming  of  this 
problem.  Therefore,  one  may  consider  this  velocity  maximum  to  be  rather  constraed.  However, 
for  simulations  of  physical  phenomena,  the  development  of  new  velocity  maxima  may  be 
misleading  or  even  catastrophic. 

The  shock  mbe  is  also  a  simple  problem,  but  it  exhibits  some  in^rtant  fluid  phenomena.  The 
domain  is  also  a  long  straight  tube,  or  chamber,  which  is  divided  by  a  diaphragm.  The  initial 
mass  density  on  one  side  of  the  diaphragm  is  larger  than  what  is  on  the  other  side,  but  both  sides 
have  the  same  initial  temperature.  For  an  ideal  gas  equation  of  state,  the  initial  pressure  is 
directly  proportional  to  the  initial  density.  When  the  diaphragm  is  released,  the  gas  with  the 


(a)  Monotonic  momentum  difference.  (b)  Altered  momentum  difference. 

Figure  7.  Advection  of  the  step  function  of  peipendicular  velocity. 
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greater  density  will  expand  into  the  lower  density  gas,  launching  a  shock  wave  ahead  of  the 
diaphragm  and  creating  a  rarefaction  wave  behind  it.  In  so  doing,  the  pressure  equilibrates 
across  the  diaphragm.  An  analytic  solution  may  be  found  for  this  problem  using  the  method  of 
characteristics  for  the  rarefaction  wave  and  the  Rankine-Hugoniot  relations  for  the  shock 
(Ref.  10). 

Figure  8  shows  the  analytic  solution  and  three  solutions  calculated  with  MACH2  at  30  ps  for  a 
shock  tube  with  a  y  =  S/3  gas  and  an  initial  density  ratio  of  four  across  the  diiq)hragm.  The 
diaphragm  is  initially  located  at  the  0.8-m  position,  and  the  initial  sound  speed  is  12.4  km/s. 
Note  that  the  results  with  the  van  Leer  algorithm  improve  the  performance  for  the  advection  of 
the  contact  surface,  which  is  the  discontinuity  at  the  released  diaphragm,  compared  with  the 
results  from  the  Godunov  method.  The  third  MACH2  curve  is  a  Lagrangian  version  of  the  same 
shock  tube.  It  maintains  a  perfect  contact  surface,  but  performs  only  slightly  better  than  the 
Eulerian  simulations  for  the  shock  and  rarefaction  waves.  Thus,  the  error  in  modeling  the  two 
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Figure  8.  Shock  tube  results. 


waves  may  be  attributed  to  the  Lagrangian  phase  and  not  the  advection  phase.  All  three 
simulations  were  mn  with  fiiUy-advanced  time  centering  for  the  implicit  Lagrangian  algorithm. 
When  the  Lagrangian  simulation  is  repeated  with  half-advanced  time  centering,  the  waves  are 
sharper,  but  the  solution  is  also  oscillatory. 


5.0  DISCUSSION 


The  van  Leer  advection  algorithm  currently  in  MACH2  has  been  used  and  upgraded  over  the 
past  2  1/2  years.  It  has  been  used  for  many  complicated  simulations  of  hydrodynamic  and  MHD 
phenomena.  The  algorithm  seems  to  be  rather  robust  and  does  not  require  special  attention  in 
most  cases.  In  a  practical  sense  it  provides  efficiency.  When  simulating  complicated 
phenomena,  one  typically  uses  only  enough  grid  resolution  to  provide  reasonable  convergence 
towards  a  solution— hopefiiUy  the  correct  solution,  to  save  on  computation  and  personal  time. 

For  a  rough  estimate,  one  needs  about  half  as  many  cells  in  each  dimension  with  the  van  Leer 
algorithm,  in  comparison  with  what  is  needed  with  the  Godunov  method,  to  achieve  the  same 
level  of  convergence  for  dynamic  simulations.  When  one  considers  that  an  increased  cell  size 
also  increases  the  allowed  time  increment,  the  savings  in  computation  time  can  be  close  to  an 
order  of  magnitude. 

The  algorithm  presented  here  should  not  be  considered  a  final  state.  If  the  user  has  ideas  for 
improvements  or  has  special  needs  for  a  particular  problem,  the  author  encourages  him  to  pursue 
them.  Writing  code  for  MACH2  is  fairly  easy  after  one  learns  the  block  structure  and  the  tool 
routines  for  setting  boundary  conditions. 
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APPENDIX 


*dktmspt 

subroutine  tn  >t(dloginmx) 

c - invoke  the  van  Leer  transport  loop  for  mass,  energy, 

c - magnetic  Held  and  momenta. 

cdir$  nolist 

include  ’common.h’ 
include  ’inputcom.h’ 
include  ’pointer.h’ 
include  ’mgcom.h’ 
cdir$  list 

pointer(  kp006,  gntrol(  0:ip2, 0:jp2) ) 
pointer(  kpOlO,  giyrol(  0:ip2, 0:jp2) ) 

c - initialize  the  volume  integrated  quantities. 

do  1001blk=  l,nblk 
call  setblk 

call  tmsinit(dlogmmx) 

100  continue 

c - first  transport  mass: 

lblk=l 
call  setblk 

call  tmslp(’other’,rol,mp) 

c - save  the  mass  gradients  for  momentum 

do  110  lblk=l,nblk 
call  setblk 

call  bkpntrs(lblk,lblk,all,c’  11, all, cell) 
call  bkcpyvf(dquanx,dquany,grxrol,gryrol) 

110  continue 

c - mass  of  material  2: 

if  (con2on)  then 
lblk=  1 
call  setblk 

call  tmslp(’other’,con24np2) 

c - if  con2on  these  gradients  are  needed 

do  120  lblk=l,nblk 
call  setblk 

call  bkpntrsOblk,lblk,aU,ceU, all, cell) 
call  bkaddar(dquanx,grxrol,grxrol) 
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call  bkaddar(dquany,gryrol,gryrol) 

120  continue 
endif 

c — -energy: 

lblk=l 
call  setblk 

caU  tmslp(’otheT’,sel,ep) 

c - ion  energy 

if  ( tsplit  .ne.  0 )  then 
lblk=  1 
call  setblk 

call  tmslp(’other’,sieion,eip) 
endif 

if  (strength)  then 
lblk=  1 
caU  setblk 

call  tmslp< ’other ’,sigdxxl,sigdxx) 

lblk=l 

call  setblk 

call  tnislp(’other’,sigdxyl,sigdxy) 

lblk-1 

call  setblk 

call  tmslp<’other’,sigdxzl,sigdxz) 
lblk=  1 
call  setblk 

call  tnislp(’othcr’,sigdyyl,sigdyy) 
lblk=  1 
call  setblk 

call  tmslp(’other’,sigdyzl,sigdyz) 
endif 

c - vertex  centered  momenta  get  special  attention  for 

c - boundary  conditions  and  fluxing. 

c - three  momentum  components: 

lblk=  1 
call  setblk 

call  tmsmm(  ul,  up ) 
lblk=  1 
call  setblk 

call  tmsmm(  vl,  vp ) 
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lblk=l 
call  setblk 

call  tmsmin(  wl,  wp ) 


c - magnetic  fields 

if  (magon)  then 

c - find  the  gradients  of  (bxl,byl,bzl) 


call  tmsbgr(brbzon) 

if  (bibzon)  then 
dol60iblk=l4iblk 
call  setblk 

c - find  the  poloidal  fields  for  fluxing 

call  tmsinib 
160  continue 
endif 

dol701blk=l^blk 


call  setblk 
if  (brbzon)  then 

c - find  E  dot  dl  for  poloidal  fluxes 

call  tmsedl(dt) 

c - transpoit  die  poloidal  flux 

call  tmsflux 

c - compute  new  poloidal  field 

call  tmsbxby 
endif 

c - transport  out-of-plane  magnetic  flux 

call  tmsbz 


170  continue 
endif 

c — -divide  returned  quantities  by  new  cell  mass  and  update  density. 
do2001blk=l4iblk 
call  setblk 
call  tmsfin 
200  continue 

return 

end 
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*dktmsinit 

subroutine  tinsinit(dlognmn) 

c - initialize  the  cell  integrals  for  advection 

cdir$  nolist 

include  ’conimon.h* 
include  ’inputcom.h’ 
include  ’pointerJi’ 
cdir$list 

dimension  dlogm(mxij) 
dimension  con2t(  0:mxij  ,  0:mxij ) 

common  /flxpnt/  istait(mxblks)4end(mzblks), 
%  jstait(mxblks)Jend(mxblks) 

ifstart  =  istart(lblk) 
ifend  =  iend(lblk) 
jfstart  =  jstart(lblk) 
jfend  =  jend(lblk) 


t3=  1./3. 
do  100  j=0  Jpl 
do  100  i=0,9l 

tnp(ij)  *ro(ij)  *  oldvol(i  j) 
mp2(ij)  =  zero 
con2t(iJ)  =  con2(iJ) 
cp(ij)  =mp(ij)  *sel(ij) 

cip(ij)  =mp(i,j)  *sieion(ij) 

sigdxx(ij)  =  mp(i,j)  *  sigdxxl(ij) 
sigdxy(i J)  =  mp(i j)  *  sigdxyl(g) 
sigdxz(i j)  =  mp(ij)  ♦  sigdxzl(ij) 
sigdyy(ij)  =  miKiJ)  *  sigdyyl(ij) 
sig(fyz(i j)  =  m^i j)  ♦  sigdyzl(i J) 
sel(i  j)  =  rol(i  j)  *  scl(i  j) 

sieion(i  j)  =  rol(ij)  *  sieion(ij) 
sigdxxl(i  j)  =  rol(i  j)  *  sigdxxl(ij) 
sigdxyl(i  j)  =  rol(ij)  ♦  sigdxyl(ij) 
sigdxzl(i  J)  =  rol(io)  *  sigdxzl(ij) 
sigdyyl(ij)  =  rol(ij)  ♦  sigdyyl(ij) 
sigdyzl(i  j)  =  rol(i  j)  ♦  sigdyzl(i  j) 
ujKij)  =0. 
vp(ij)  =0. 
wp(ij)  =0. 

c - define  a  lagrangian  volume 
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lagvol(i  j)  =  mp(i  j)  /  ( rol(ij)  +  tiny  ) 

c - volume  flux  out  of  the  bottom  of  this  cell. 

rub  =  ( 2.*r<ij)  +  r(i+l  J) )  *  url(ij)  + 

%  ( r(i  j)  +  2.*r(i+l  j) )  ♦  url(i+l  j) 

rvb  =  ( 2.*r(ij)  +  r(i+l,j) )  *  vrl(ij)  + 

%  ( r(io)  +  2.*r(i+lo) )  ♦  vrl(i+l.j) 

dxbott(i,j)  =  -0.5  *  dt  *  t3  * 

%  (  rub*(y(i+l J)  -  y(ij))  -  rvb*(x(i+l  j)  -  x(ij))  ) 

c - volume  flux  out  of  the  left  of  this  cell. 

ml  =  (  2.*r(iJ)  +  ftij+l) )  *  url(ij)  + 

%  (r(ij)  +  2.*r(ij+l))*url(ij+l) 

rvl  =  (  2.*r(ij)  +  r(ij+l) )  *  vrl(ij)  + 

%  ( r(i,j)  +  2.*r(iJ+l) )  ♦  vrl(ij+l) 

dxleft(ij)  =  -0.5  *  dt  *  t3  ♦ 

%  ( rul*(y(io)  -  y(ij+l))  -  rvl*(x(ij)  -  x(ij+l)) ) 

100  continue 

c - separate  con2  from  coni  for  fluxing  to  keep  con2  <  1. 

if  (con2on)  then 
do  200  i=0,4)l 
do  200  j=0,jpl 

mp2(i,j)  =  con2t(i,j)  *  mp(ij) 
mp(i,j)  *  ( 1-  -  con2t(ij) )  *  mp(i,j) 
con2(iJ)  =  con2t(ij)  ♦  rol(ij) 
rol(i,j)  =  ( 1.  -  con2t(ij) )  ♦  rol(ij) 

200  continue 
endif 

c - time  step  control  based  on  volume  flux  to  cell  volume  ratio— compare 

c - with  both  new  volumes  and  lagrangian  volumes. 

do  300  j=jfstart  Jfend 
do  350  i=l  ,icels 

volmin  =  min(  lagvol(i,j),  lagvol(ij-l), 

%  one  /  rvol(ij),  1.  /  rvol(ij-l) ) 

dlogm(i)  =  dxbott(ij)  /  volmin 
350  continue 

if  (j  .ne.  jpl)  then 
idlogm  =  isamax(icels,dlogm(l),l) 
if  ( abs(  dlogm(idlogm) )  .gt.  dlogmmx  )  then 
idtc  =  idlogm 
jdtc  =  -j 
Idtc  =  Iblk 

dlogmmx  =  abs(  dlogm(idlogm) ) 
endif 
endif 
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300  continue 


do  400  i=ifstart^end 
do450  j=l  jccls 

volmin  =  min(  lagvol(ij),  lagvol(i-l  j), 

%  one  /  rvol(i  J),  one  /  rvol(i-l  j) ) 

dlogm(j)  =  dxleft(ij)  /  volmin 
450  continue 

if  (i  .ne.  ^1)  then 
jdlogm  =  isamax(jcels,dlogm(l),l) 
if  ( abs(  dlogm(jdlogm) )  .gt.  dlogmmx )  then 
jdtc  =  jdlogm 
idtc  =  -i 
Idtc  =  Iblk 

dlogmmx  =  abs(  dlogm(jdlogm) ) 
endif 
endif 

400  continue 

return 

end 
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*dktinslp 

subroutine  tmslp(quantype,tnisfTom,tmsto) 

c - transport  the  quantity  tmsfrom 

c - using  the  van  Leer  transport  scheme  where  the  flux 

c - across  boundary  i+1  for  a  positive  velocity,  v,  is 

c - v(  q(i)  +  .5(1  -  V  dt/dx)  X  dq(i) ) ,  and 

c - dq(i)  =  ( q(i+l)  -  q(i-l) )  /  2  . 

c - this  first  loops  over  the  blocks  to  find  the  dq’s, 

c - then  loops  over  the  blocks  to  transport  q. 

cdir$  nolist 

include  ’patamcom.h* 
include  ’inputcom.h’ 
include  ’pointerJi’ 
cdir$  list 

character*(*)  quantype 

dimension  tms£rom(0:^2,0:jp2),  tmsto(0:ip2,0:jp2) 
pointer  ( kpptl,  qul(0:q)2, 0:jp2) ) 
pointer  ( Iq^m,  qun(0:ip2, 0;jp2) ) 

c - get  pointer  numbers  for  input  arrays 

ikpptl  =  lindex  (tmsfrom) 
ilq)pm  =s  lindex  (tmsto) 

c - gradient  loop 

if  (quantype  ne.  ’momentum’)  then 

do  100  lblk=l,nblk 
call  setblk 

kpptl  =  lpoint(  ikpptl,  Iblk) 
call  tmsgr(quantype,qul) 

100  continue 

c - transport  loop 

do  200  lblk=l,nblk 
call  setblk 

kpptl  =  lpoint(  ikpptl,  Iblk) 
kpptn  =:  lpoint(  ikpptn,  Iblk) 
call  tmsdqbc 
call  tmsadv(qul,qun) 

200  continue 
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else 


do  300  lblk=l^blk 
call  setblk 

kpptl  =  lpoint(  ilq^tl,  Iblk) 
kpptn  =  lpoint(  ikpptn,  Iblk) 
call  tnismmgr(qul) 
call  tinsadv(qul,qun) 

300  continue 
endif 

return 

end 
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*dktnisgr 

subroutine  tmsgr(quantype,quan) 

c - find  the  limited,  centered  differences  for  use  in  the 

c - van  Leer  transport  scheme. 

cdir$  nolist 

include  ’common.h’ 
include  ’inputcom.h’ 
include  ’pointer.h’ 
cdir$  list 

character*(*)  quantype 

common  /flxpnt/  istait(mxblks),iend(mxblks), 

%  jstait(mxblks)Jend(mxblks) 

dimension  quan(0:q)2,0:jp2) 

do  100  j=0  jp2 
do  100  i=0,ip2 
dquanx(ij)  =  0. 
dquany(ij)  =  0. 

100  continue 

c - set  gradient  ranges;  limit  it  along  walls  to 

c - prevent  confusion  (taken  care  of  in  rmshbcs), 

c - forcing  donor  cell  there  except  poloidal  B. 

if  ( .not.  donor(lblk) )  then 
if  (quantype  .eq.  ’polbfld’)  then 
igrxst  =  1 
igrxend  >  icels 
jgrxst  =  0 
jgrxend  =  jpl 
igiyst  =  0 
igryend  =  ipl 
jgryst=  1 
jgryend  =  jcels 
else 

igrxst  =  istart(lblk) 
igrxend  =  iend(lblk)  - 1 
jgrxst=  1 
jgrxend  =  jcels 
igryst=  1 
igryend  =  icels 
jgryst  =  jstait(lblk) 
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jgryend  =  jend(lblk)  - 1 
endtf 

do  200  j=jgiyst  jgryend 
do  200  i=igTyst4giyend 
diffb  =  quan(ij)  -  quan(ij-l) 
difft  =  quan(ij-f  1)  -  quan(ij) 
diffc  =  ( quan(ij+l)  -  quan(i  j-1) )  /  2.d0 
sdiffb  =  sign(  one,  diffb ) 
sdifft  =  sign(  one,  difft ) 
sdiffc  =  sign(  one,  diffc  ) 
dlimb  =  diffb  *  2.d0 
dlimt  =  difft  *  2.d0 

sdqy  =  max(  zero,  sdiffb  *  sdifft )  /  sdiffc 

dquany(i  j)  =  sdqy  ♦  min(abs(dlimb),abs(dUrnt),abs(diffc)) 

srmrvl  =  sign(  one  ,  ( ro(i  j)  -  rofvl ) ) 

grmlt  =  max(  zero ,  srmrvl ) 

dquany(i  j)  =  grmlt  ♦  dquany(i  j) 

200  continue 

do  300  i=igrxst,igrxend 
do  300  j=jgrxst  jgrxend 
diffl  =  quan(i,j)  -  quan(i-l  j) 
difft  =  quan(i+l  j)  -  quan(ij) 
diffc  =  ( quan(i+l  j)  -  quan(i-l  j) )  /  2.d0 
sdiffl  s  sign(  one,  diffl ) 
sdifft  =  sign(  one,  difft ) 
sdiffc  =  sign(  one,  diffc  ) 
dliml  =  diffl  *  2.d0 
dlimr  =  difft  *  2.d0 

sdqx  =  max(  zero,  sdiffl  *  sdifft )  /  sdiffc 

dquanx(i  j)  =  sdqx  ♦  min(abs(dliml),abs(dlimr),abs(diffc)) 

srmrvl  =  sign(  one  ,  ( ro(i  j)  -  rofvl ) ) 

grmlt  =  max(  zero ,  srmrvl ) 

dquanx(i  j)  =  grmlt  *  dquanx(i  j) 

300  continue 
endif 

return 

end 
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*dktmsdqbc 

subroutine  tmsdqbc 

c - communicates  the  van  Leer  gradient  across 

c - block  boundaries. 

cdir$  nolist 

include  ’common.h’ 
include  ’pointer .h’ 
include  ’bccommon.h’ 
include  ’inputcom.h’ 
cdir$list 


do  100  i=l,4 
ibdry  =  iproseq(Ublk) 

Inbr  =  knbr(ibdiy4blk) 
if  (Inbr  .ne.  0)  then 
call  sembrb(lnbr) 

c  from  to  range 

call  bcpntrs(ibdiy,nebr, edge, this, ghst.edge, cell) 
call  bccpyvf(dqxnbr,dqynbt,dqulx,dquly) 
endif 

100  continue 

return 

end 
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*dktmsadv 

subroutine  tmsadv(quanf, quant) 

c - transport  a  quantity  with  flux  based  on  what  is  in 

c - quanf  into  what  is  in  quant  with  a  van  Leer  scheme. 

c - be  aware  that  what  come  into  this  routine  is  something 

c - per  unit  volume  (something  X  density)  and  it  returns  a 

c - volume  integrated  quantity  (to  be  divided  by  new  cell  mass). 

cdir$  nolist 

include  ’common.h’ 
include  ’inputcom.h’ 
include  ’pointer.h’ 
cdir$  list 

dimension  quanf(0:ip2, 0:jp2),  quant(0:ip2, 0:jp2) 

common  /flbipnt/  istart(mxblks)4end(mxblks), 

%  jstart(mxblks)jend(mxblks) 

ikpptt  =  lindex(quant) 
ifstart  =  istart(lblk) 
ifend  =  iend(lblk) 
jfstart  =  jstait(lblk) 
jfend  =  jend(lblk) 


c - flux  in  j-direction 


do  100  j=jfstart  Jfend 

c - compute  the  fluxes  between  this  row  and  the  row  below. 

do  ISO  i=l,icels 

c - advection  (don’t  forget  dxbott  >  0  implies  downward  flow) 

sdv  =  sign(  one  ,  dxbott(i,j) ) 

sigma  =  O.SdO  ♦  ( (one  +  sdv)  ♦  dxbott(iJ)/lagvol(i  j) 

%  +  (one  -  sdv)  ♦  dxbott(ij)Aagvol(i,j-l) ) 

qubsp  =  ( quanf(ij)  -  O.SdO  *  (one  -  sigma)  *  dquany(ij) ) 
qubsm  =  ( quanf(i,j-l)  +  O.SdO  *  (one+sigma)  *  dquany(ij-l) ) 
dqbs  =  O.SdO  *  dxbott(i,j)  ♦ 

%  ( (one  +sdv)  ♦  qubsp  +  (one  -  sdv)  ♦  qubsm  ) 

quant(ij)  =  quant(ij)  -  dqbs 
quart(ij-l)  =  quant(ij-l)  +  dqbs 

ISO  continue 
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100  continue 


c - flux  in  i-direction 


do  200  i=ifstart4fend 

c - compute  the  fluxes  between  this  column  and  the  colunm  to  the  left. 

do  250  j=l  Jcels 

c - advection  (don’t  forget  dxleft  >  0  implies  flow  to  ihe  left) 

sdv  =  sign(  one  ,  dxleft(i  J) ) 

sigma  =  0.5d0  *  ( (one  +  sdv)  *  dxleft(i j)Aagvol(iJ) 

%  +  (one  -  sdv)  *  dxleft(ij)Aagvol(i-l  j) ) 

qulsp  =  ( quanf(ij)  -  0.5d0  ♦  (one  -  sigma)  *  dquanx(ij) ) 
qulsm  =  ( quanf(i-l  J)  +  0.5d0  *  (one+sigma)  ♦  dquanx(i-l  J) ) 
dqls  =  0.5d0  *  dxleft(ij)  * 

%  ( (one  +  sdv)  •  q^sp  +  (one  -  sdv)  *  qulsm ) 

quant(ij)  =  quant(ij)  -  dqls 
quant(i>l  J)  =  quant(i-l  J)  +  dqls 

250  continue 

200  continue 

return 

end 
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*dk  tmsinin 

subroutine  tmsnim(tmsfrom,tinsto) 

c - Use  the  mach2  momentum  con^nents  [cell-ro  *  vert-v] 

c - and  use  only  v  *  grad(ro) 

c - to  avoid  advecting  a  lot  of  mass  and  little  momentum. 

cdir$  nolist 

include  ’paramcom.h’ 
include  ’pointer.h’ 
cdir$  list 

dimension  tmsfrom(0:^2,0:jp2),  tmsto(0:ip2,0:jp2) 
pointer(  kpptf,  qul(0:ip2, 0:jp2) ) 
pointer(  kpptt,  qun(0:ip2, 0:jp2) ) 

c - get  pointer  numbers 

ikpptf  =  lindex(tmsffom) 
ikpptt  =  lindex(tmsto) 

do  1001blk=  l^blk 
call  setblk 

kpptf  =  lpoint(  ikpptf,  Iblk ) 
call  tmsmmin(qul) 

100  continue 

c - advect  the  four  ’moments’: 

Iblk  =  1 
call  setblk 

call  tmslp(’momcntum’40wl,rowlt) 
lblk=  1 
call  setblk 

call  tmslp(’momentiun’,row24’ow2t) 
lblk=  1 
call  setblk 

call  tmslp(’momentum’,row3,row3t) 

Iblk  =  1 
call  setblk 

call  tmslp{’momentum’4X)w4,row4t) 

c - recreate  vertex  momenta. 

do2001blk=  l,nblk 
call  setblk 

kpptt  =  lpoint(  ikpptt,  Iblk  ) 
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call  tinsininf(qiin4]q)ptt) 
200  continue 

return 

end 


^dktmsmmin 

subroutine  tmsnimin(vel) 

c - initialize  the  four  ’moments’  of  momentum  for  one  component. 

cdir$  nolist 

include  ’common.h’ 
include  ’inputcom.h’ 
include  ’pointer.h’ 
cdir$list 

dimension  vel(0:ip2, 0:jp2) 

do  100  j=0,jpl 
do  100  i=0,q)l 
vml  =  vel(i+l  j) 
vm2  =  vel(i+l  j+1) 
vm3  =  vel(ij+l) 
vm4  =  vcl(i  j) 
tmass  =  ro(i  J)  *  oldvol(i  J) 
rowl(i,j)  =  vml 
row2(i,j)  =  vm2 
row3(i  j)  =  vm3 
row4(i  j)  =  vni4 
rowlt(i  j)  s  tmass  *  vml 
rovv2t(i  j)  =  tmass  *  vm2 
rovv3t(i  j)  =:  tmass  ♦  vm3 
rovv4t(i  j)  =  tmass  *  vm4 
100  continue 

return 

end 
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*dktinsininf 

subroutine  tmsnunf(velp4kpptt) 

cdir$  nolist 

include  ’conimon.h’ 
include  ’inputconi.h’ 
include  ’pointerJi’ 
cdiT$list 

dimension  velp(0:ip2, 0:jp2) 
pointer(npptt,  velpnbr(0:inp2, 0:jnp2) ) 


c - acquire  the  part  of  the  vertex  momentum  along  boundaries 

c - that  was  created  in  previous  blocks. 


do  100  ibdry=l,4 
Inbr  =  knbifibdryJLblk) 
if  (Iblk  .gt.  Inbr  .and.  Inbr  .gt.  0)  then 
call  setnbrb(lnbr) 
npptt  =  lpoint(  ikpptt,  Inbr ) 
c  from  to  range 

call  bcpntrs(ibdiy^ebr,edge,this,edge,edge,veit) 
call  bccpysc(velpnbr,vclp) 

elseif  Gblk  .eq.  Inbr  .and.  ibdry  .gt.  nbibdy(ibdiy4blk))  dien 
call  sembrb(lnbr) 
npptt  =  lpoint(  ikpptt,  Inbr ) 
c  from  to  range 

call  bcpntrs(ibdry,nebr,edge,this, edge, edge, vert) 
call  bccpysc(velpnbr,velp) 
endif 

100  continue 

do  200  icmr=l,4 
Idnbr  =  ldignbr(icmr4blk) 
if  (Iblk  .gt.  Idnbr  .and.  Idnbr  .gt.  0)  then 
call  sembrb(ldnbr) 
npptt  =  lpoint(  ikpptt,  Idnbr ) 
c  from  to 

call  ccpntrs(icmr4iebr,edge,0,this,edge,0,vert) 
call  cccpysc(velpnbr,velp) 
endif 

200  continue 

c - recombine  the  vertex  momentum  m  this  block. 

do  375  j=l  jcels 


38 


do  300  i=l,icels 

velp(i+l  J)  =  velp(i+l J)  +  0.25  ♦  rowlt(ij) 

300  continue 
do  325  i=l,iccls 

velp(i+l  j+1)  =  velp(i+l  J+1)  +  0.25  *  row2t(iJ) 

325  continue 
do  350  i=l,icels 

velp(ij+l)  =  velp(ij+l)  +  0.25  *  iDw3t(ij) 

350  continue 
do  375  i=l,icels 

velp(ij)  =  velp(ij)  +  0.25  *  row4t(i j) 

375  continue 

c - ^put  the  contribution  of  this  block  into  the  boundaries 

c — -of  previous  blocks. 

do  400  ibdry=l,4 
Inbr  =  knbr(ibdrydblk) 
if  (Iblk  .gt.  Inbr  .and.  inbr  .gt.  0)  then 
call  setnbrb(lnbr) 
npptt  =  lpoint(  ikpptt,  Inbr ) 
c  from  to  range 

call  bcpntrs(ibdry, this, edge^ebr, edge, edge,vert) 
call  bccpysc(velp,velpnbr) 

elseif  Gblk  .eq.  Inbr  .and.  ibdry  .gt.  nbtbdy(ibdry4blk))  then 
call  setnbrb(lnbr) 
npptt  =  lpoint(  ikpptt,  Inbr ) 
c  frx)m  to  range 

call  bcpntrs(ibdry,dus,edge4iebr,edge,edge,vert) 
call  bccpysc(velp,velpnbr) 
endif 

400  continue 

do  500  icnir=l,4 
Idnbr  =  Idignbrficmrdblk) 
if  (Iblk  .gt.  Idnbr  .and.  Idnbr  .gt.  0)  then 
call  setnbrb(ldnbr) 
npptt  =  lpoint(  ikpptt,  Idnbr ) 
c  from  to 

call  ccpntrs(icmr,this,edge,0,nebr,edge,0,vert) 
call  cccpysc(velp,velpnbr) 
endif 

5(K)  continue 

return 

end 
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